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Abstract 

In this paper, we propose a novel method for parcellating the human brain into 193 anatomical structures based on 
diffusion tensor images (DTIs). This was accomplished in the setting of multi-contrast diffeomorphic likelihood fusion using 
multiple DTI atlases. DTI images are modeled as high dimensional fields, with each voxel exhibiting a vector valued feature 
comprising of mean diffusivity (MD), fractional anisotropy (FA), and fiber angle. For each structure, the probability 
distribution of each element in the feature vector is modeled as a mixture of Gaussians, the parameters of which are 
estimated from the labeled atlases. The structure-specific feature vector is then used to parcellate the test image. For each 
atlas, a likelihood is iteratively computed based on the structure-specific vector feature. The likelihoods from multiple 
atlases are then fused. The updating and fusing of the likelihoods is achieved based on the expectation-maximization (EM) 
algorithm for maximum a posteriori (IVIAP) estimation problems. We first demonstrate the performance of the algorithm by 
examining the parcellation accuracy of 18 structures from 25 subjects with a varying degree of structural abnormality. Dice 
values ranging 0.8-0.9 were obtained. In addition, strong correlation was found between the volume size of the automated 
and the manual parcellation. Then, we present scan-rescan reproducibility based on another dataset of 16 DTI images - an 
average of 3.73%, 1.91%, and 1.79% for volume, mean FA, and mean MD respectively. Finally, the range of anatomical 
variability in the normal population was quantified for each structure. 
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introduction 

For quantitative analysis of the liuman brain anatomy, defining 
structures or regions of interest (ROIs) is one of the first essential 
steps. There are many types of automated or manual approaches 
that have been proposed to define ROIs in the brain, based on 
locations and contrasts of the structures. These methods often 
require a priori knowledge as a form of atlas. For manual ROl 
drawing, an atias could be a simple pictorial representation of a 
structure of interest, which guides operators to define the 
boundary. The manual delineation, while often used as a gold 
standard, is a time-consuming approach. Various types of 
automated parcellation tools have been proposed, which try to 
define the boundary of anatomical structures based upon image 
contrasts [1-17]. Some of the advanced tools incorporate a. priori 
knowledge about the location of the target structures as a form of 
probabilistic atlas [18-26]. This location constraint prevents the 
contrast-based boundary definition from leaking into unlikely 
regions. 

To use the probabilistic location information in the atlas, the 
adas has to be registered, or warped, to each subject image, in 
which voxel-to-voxel correspondence is established between the 
two coordinate systems where the atlas image and the subject 



image are defined. The concept of brain mapping also leads to an 
alternative approach for automated structural parcellation, which 
is called adas-based parcellation [3,7,11,27-35]. Namely, if the 
voxel-to-voxel mapping is perfecdy accurate, any arbitrary 
structures can be defined only once in the adas and such 
anatomical definitions can be transferred to the images being 
mapped to. There are no constraints in the number of structures 
or the way the structures are defined in the atlas. In this sense, we 
can assume that the whole-brain voxel mapping inherendy 
contains parcellation tools for potentially all definable structures 
inside the brain. 

The adas-based parcellation, however, is accurate only if the 
voxel-to-voxel brain mapping correctly defines corresponding 
voxels between the two images, which is not always guaranteed. 
The accuracy level is influenced by various sources of differences 
between the atias and the subject brains; these could be 
morphological (atrophy, hypertrophy, malformation, etc.) or 
contrast (biological such as signal hyperintensity or hypointensity 
or procedural such as imaging parameters). 

To reduce the impact of erroneous mapping of voxels, and 
consequent mis-parceUation of target structures, multi-atlas 
approaches have been postulated. Suppose that the hippocampus 
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is defined in JVdifTerent atlases and each atlas is warped to a to-be- 
parcellated subject image, then the A^definitions (labels) of the 
hippocampus are casted to the subject space, which can be fused 
(a.k.a "label fusion") based upon pre-defined algorithms such as 
those proposed by [7,36-42]. If the mis-registration of each atlas 
causes random errors, the errors should be reduced by integrating 
A'^definitions. It has been shown that simple label fusion techniques 
based on majority voting yield robust parcellations [7,41,43]. 
More recendy, weighted majority voting strategies, by incorpo- 
rating intensity information, demonstrated significant improve- 
ment in the parcellation accuracy. A variety of weighting 
approaches, based on intensity similarity metrics, have been 
proposed - global [36], local [44,45], semi-local [45,46], and non- 
local [47]. In addition to voting, a statistical fusion technique (i.e. 
Simultaneous Truth and Performance Level Estimation, STAPLE 
[42]) and a collection of its variants [48-50] have been proposed, 
in which a stochastic model of rater behavior has been 
incorporated in the estimation process. Compared with voting 
techniques, the main limitation of statistical fusion strategies is that 
the decision rule is independent of the image intensity while the 
major advantage is its underlying elegant mathematical theory. 
Initial attempts to incorporate the intensity information into the 
STAPLE framework rely on a priori similarity measures [51,52] or 
estimating the voxelwise correspondence between the registered 
rater and the subject using intensity information [49]. 

We have recently introduced the diffeomorphic likelihood 
fusion algorithm (DLFA) as an approach to integrate anatomical 
information from multiple T 1 -^\(ught(xi atlases and fuse their 
anatomical features [53,54]. Unlike previous label fusion algo- 
rithms, DLFA does not fuse a set of binary label maps obtained 
from the atias-to-subject propagations. DLFA poses the parcella- 
tion problem in the framework of maximum a posteriori (MAP) 
estimation, estimating the maximizing parcellation labels given the 
observable image intensity, similar to the idea proposed in [6]. 
The MAP estimation is handled within the class of generative 
models by representing the observable imagery as a conditionally 
Gaussian mixture random field, conditioned on the random atlas- 
label pair and the diffeomorphic change of coordinates for each 
label. The atias-label pair and their diffeomorphic correspondenc- 
es are unknown and viewed as latent variables. Locality is 
introduced into the global representations of the deformable 
templates by allowing different atlas-label pairs to be used to 
interpret different voxels or different structures, under the 
assumption that the local optimal diffeomorphism varies from 
label to label for a given atias. The MAP estimation is solved by 
iterating between fixing the local optimal diffeomorphisms and 
obtaining the maximizing parcellation labels, and then locally 
optimizing the local diffeomorphisms for the fixed parcellation, in 
an EM fashion [55]. The atias-dependent structure-specific local 
diffeomorphisms are estimated in the E-step in the EM algorithm. 

The purpose of this paper is to extend the DLFA to diffusion 
tensor imaging by incorporating multiple-contrast information. It 
arises naturally to extend a single contrast image such as Tl- 
weighted images to multi-contrast images (e.g. eigenvalues and 
eigenvectors of DTI) by assuming conditional independence in 
computing p{l\ W,a), where Idenotes the measurable image, W 
denotes a given parcellation label, and a the randomly selected 
atlas-label pair. Previously, vector-to-vector or tensor-to-tensor 
registration algorithms have also been introduced [56-59], which 
was further extended to multi-channel image registration, in which 
multiple-contrast information, such as FA, diffusivity, and fiber 
orientation, is used simultaneously to drive the registration 
algorithm [60]. These ideas could improve the registration 
accuracy between each atias and the testing subject, but the 



incorporation of the multiple-contrast information in the multi- 
atias likelihood fusion process has not been introduced so far. 

In this paper, we introduce a framework to incorporate the 
multi-contrast intensity information generated in DTI into the 
multi-atlas DLFA framework and apply it to whole brain 
parcellations into 159 structures. In the Tl case, the distribution 
of the intensity in each structure is modeled as a single Gaussian. 
In the DTI case, we use five intensity elements ([FA, MD, and 
fiber angle (a unit vector)]) with the intensity distribution of each, 
in every single label, being modeled as a Gaussian Mixture Model 
(GMM), the parameters of which are computed using maximum- 
likelihood estimation. In this study, we examine the parcellation 
accuracy of the method on 25 patient data with a varying degree 
of pathology. We also present the scan-rescan reproducibility of 
the method on another dataset of 16 healthy subjects which were 
scanned twice. In addition, the ranges of anatomical variability of 
all the structures in the 16 subjects were characterized. 

Methods and Materials 

Patient populations 

AU subjects used as atiases and the first testing dataset were 
obtained from the existing clinical database of pediatric brain 
MRI, and were older than 24 months of age. DTIs from sixteen 
subjects (Female = 7, Male = 9, age = 7.67+/ — 4.12) were used to 
create the multiple atiases (Table 1). Among these sixteen subjects, 
ten subjects were diagnosed as normal. In order to cover the wide 
range of anatomical phenotypes in the multiple atiases, 6 cases 
with different types of anatomical abnormalities were also included 
in the atlas set as shown in Table 1. 

The accuracy of the parcellation obtained from multi-atias 
DLFA was tested using 25 patients (Female =10, Male =15, 
age = 7.88+/-4.80). As tabulated in Table 2, 10 subjects (Test 
#1-#10) presented a normal MR anatomy and the other 15 
subjects pr(;scnted a variety of anatomical abnormalities; seven 
(Test #11-#17) were evaluated as mild to moderate anatomical 
change and eight (Test #18-#25) as severe abnormality based on 
a pediatric neuroradiologist's (S.Y.) visual evaluation. 

For the scan-rescan reproducibility test, sixteen healthy 
volunteers with no history of neurological conditions (8 M/8 F, 
22-61 years old, mean: 31 years old) participated in this study. 
This is the same data used by [61], where details of the protocol 
can be found. 

Ethics Statement 

This study was approved by the Johns Hopkins Medicine 
Institutional Review Board (JHM-IRB). All subjects provided 
written, informed consent for participation in accordance with the 
oversight of the JHM-IRB. 

MRI scans 

For the first dataset, MR imaging was performed using a 1 .5T 
scanner (Avanto; Siemens, Erlagen, Germany). All patients 
underwent routine clinical multiplanar Tl, T2, and FLAIR pulse 
sequences, including DTI. The DTI was obtained using a single- 
shot EPI with parallel acquisition. Diffusion weighting was 
performed along 21 independent axes with b= 1000 s/mm2, 
and repeated twice to enhance the SNR (TE = 84 ms, 
TR = 7700 ms). DTI was scanned in the axial orientation with 
an imaging matrix of 96x96 (to 192x192 with zero-filled 
interpolation), FOV 240x240, and slice thickness 2.5 mm. 

For the second dataset, subjects were scanned twice using a 3T 
MR scanner (Achieva, Philips Healthcare, Best, The Netherlands). 
The DTI dataset was acquired using a multi-slice, single-shot. 
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Table 1. 


Anatomical changes in the sixteen subjects used as the multiple atlases. 






Atlas ID 


Radiological findings 


Radiological diagnosis 


1 


White matter T2 hyperintensity involving the bilateral periventricular and deep white 
matter with restricted diffusion spots 


Drug-induced leukoencephalopathy 


2 


T2-hyperintense lesions in periventricular and subcortical white matter 


Multiple sclerosis 


3 


Multiple encephalomalacia/gliosis change related to sequela from prior ischemic events 


Moyamoya-disease 


4 


Diffuse CSF space dilatation 


Associated finding with achondroplasia 


5 


Multiple T2-hyperintense lesions in white matter and gray matter 


Neurofibromatosis typel 


6 


Mild ventricular dilatation with irregular shape and volume loss of periventricular white 
matter with posterior dominant 


Periventricular leukomalacia 


7-16 


No abnormal finding 


Diagnosed as normal 


doi:1 0.1 371 /journal.pone.0096985.t001 



echo-planar imaging (EPI), spin-echo sequence (TR/TE = 6281/ 
67 ms, SENSE factor = 2.5). Sixty-five transverse slices were 
acquired parallel to the line connecting the anterior commissure 
(AC) to the posterior commissure (PC) with no slice gap and 
2.2 mm nominal isotropic resolution (FOV = 212x212, data 
matrix = 96x96, reconstructed to 256x256). 

DTI processing 

All DTI datasets were processed offline using DTIStudio 
software (H. Jiang and S. Mori, Johns Hopkins University, 
Kennedy Krieger Institute, lbam.med.jhmi.edu or www. 
MriStudio.org) [62]. The raw diffusion-weighted images were fu-st 
co-registered to one of the bO images with a 12-parameter affme 
transformation using Automated Image Registration (AIR) [63]. 
The six elements of the diffusion tensor, the fractional anisotropy 
(FA), and the mean diffusivity (MD) were calculated. 



Initial creation of multiple atlases 

For the sixteen subjects that were selected to be the multiple 
atlases, the images were first normalized to MNI coordinates with 
a nine-parameter affine transformation. The initial parcellation of 
the brain into 159 structures was performed using the atlas-based 
automated image parcellation pipeline as described in our 
previous pubUcation [60]. We used our single-subject Eve atlas 
[64] and the accompanied brain parcellation map with 159 
structural definitions as the template, which was warped to the 16 
subjects using the three-contrast large deformation diffeomorphic 
metric mapping (LDDMM) [64—66]. The three contrasts included 
FA, MD, and the manually-delineated lateral ventricles. In a 
previous study, we tested the accuracy of this automated structural 
parcellation approach in cerebral palsy patients and excellent 
accuracy was reported [65]. In this study, we included patients 
with more severe abnormalities. If gross parcellation errors 
occurred, they were manually corrected to establish the multiple 



Table 2. Anatomical changes in 25 subjects for testing the multiple atlases application. 




Test Subject ID 


Radiological findings 


Radiological diagnosis 


1-10 


No abnormal finding 


Diagnosed as normal 


11 


Mild deep white matter T2-hyperintense change and ventricle enlargement 


Adrenoleukodystrophy 


12 


Right hemiatrophy, ventricle dilatation and mild T2-hyperintense change in deep 
gray matter 


Chronic ischemic insult 


13 


Diffuse CSF space and ventricle dilatation 


Associated finding with achondroplasia 


14 


Mild ventricle dilatation 


Associated finding with achondroplasia 


15 


T2-hyperintense lesions in periventricular and subcortical white matter and mild 
ventricle enlargement 


Multiple sclerosis 


16 


Porencephalic left ventricle dilatation and volume loss of left corticospinal tract 


Prenatal hemorrhagic insult 


17 


Asymmetrical ventricle dilatation (right>left) 


Associated finding with achondroplasia 


18 


Ventricle enlargement with multiple T2-hyperintense lesions in white matter 


Multiple sclerosis 


19 


Ventriculomegaly associated periventricular volume loss (right>left) and T2- 
hyperintense change 


Perinatal hypoischemic injury 


20 


Lateral ventricular enlargement with periventricular white matter volume loss 
and T2-hyperintense change 


Periventricular leukomalacia 


21 


Right ventricle enlargement associated with periventricular white matter volume loss 


Prenatal intraventricular hemorrhage 


22 


Diffuse parenchymal volume loss, CSF space dilatation and multiple ischemic lesions 


Congenital metabolic disease 


23 


Ventriculomegaly and thinning of corpus callosm 


Ventriculomegaly 


24 


Left hemiatrophy 


Sturge- Weber syndrome 


25 


Left parenchyamal volume loss, gliosis and lateral ventricle enlargement 


Perinatal stroke 


doi:1 0.1 371 /journal.pone.0096985.t002 
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atlases with accurate structural definitions. The non-linear image 
transformation and the atlas-based parcellation were performed 
using DifiFeoMap and RoiEditor (http://www.MriStudio.org, 

Kennedy Krieger Institute and Johns Hopkins University, X. Li, 
H. Jiang, and S. Mori). The atlas data are available at http:// 
lbam.med.jhu.edu. 

Multi-contrast likelihood-fusion 

Let {(r\W'"),(r\W"^),...Xl'"',W"'')}denotc^ N DTI atlas- 
label pairs, where N=16 in this study. Instead of using single- 
valued (T 1 -weighted) images, we use vector-valued images for 
both the atlases and the test subjects. For each atlas-label pair a, 
l'' = [lFA^lMi)J^'Iy'IzUe{aua2,-,au}, where I^j^dcnotes the 
gray-scale FA image of the atlas-label pair, denotes the gray- 
scale MD image of the atias-label pair, and I", 1^,1" denote the 
absolute values of the three elements of the primary eigenvector. 
In this sense, the image intensit)' at each voxel is a 5-element 
vector \"{x) : xeQ.-*U^ , with £1 <= being a finite grid where the 
images are defined. For the label image W in each atias-label 
pair, we define it as a function from the image domain fi to a 
subset of the non-negative integers fF''(x) : X6Q^{0,1,2,3,..., 
159}, where W''(x) = Q for voxel x belonging to the unlabeled 
background, and W''ix) = k,ke{\,2,'i,...,\59} for voxel x lal)eled 
as the k-th structure such as the left caudate, the right putamen, 
and so on. Correspondingly, we denote the to-be-parceUated test 
subject as (I, W), where I = [IvKjMaJyiJyJz\ and W is the label 
image we aim to obtain. 

For multi-contrast, multi-atias parcellation, the goal is to 
estimate the label map W associated with the im- 
agel= [/pA^^MD A^■'y.^z] of the test subject, for which we solve 
via the Maximum a Posteriori (MAP) estimation 



W= argmax/>(f^|I) = arg max/'( IF,I) . 
w w 



(1) 



fr"™=argmax Q{W;W'''^) 



(5) 



is monotonic in the incomplete data likelihood with atlas selector 
PA(x){a\l,W°''^), the proof of which can be found in [53]. 
The algorithm can be summarized as: 

Stepl: Initialize the dififeomorphism for each voxel x to be 
identical everywhere, as: ip"l'' (x) = (p^^\ Initialize IF"'''. 
Step2: Compute the approximated atias-label selector as: 

Pa(.MI{x),W'''\x)) 

_ pia,ff(x)\I(x),W<"^(x)) 



EpMri^MxlWo'dix)) 

a 

piI(x),W""(x)\a:r"(x))n(4>f(x)\a)n(a) 
EpiIix),W''Hx)\a,fJ\x)M€"i^)\«M'') 



(6) 



Step3: Obtain a new parcellation image for the test image via 
j^new = argmax QiW; W"), where QiW; W''''')k computed as: 



Q(W; W<''')=J2Y^PAid«\KxW'''''ix))\ogpilixW(x)\a) .(7) 

xeSl a 

Step4: Recalculate the diffeomorphisms of the atiases onto the 
parcellation labels via: 



q)"J''' = a.rgmaxp(l(x),W""'\a,g>)n((p\a) 
<p 

= aTgmaxp{l(x)\W"",a,q))p{W"'"'\a,^)ni(p\a) 



(8) 



To achieve this goal, we use the EM algorithm by introducing the 
latent variable ^e{ai,a2,...,ajv}that designates the random atias- 
label pair. The Qjfunction in the EM algorithm computes the log- 
likelihood of the complete data log/»(I, W,A) given the incomplete 
data — the to-be-parceUated measured image I and the previous 
parceUation label W''', 



Q(W- W"''') = E^^^^^^M){\ogP^, W\A)\\, W"} , (2) 



where 



p(l, W\A)=IIII p{lix), W(x)\af^(^^^''^ (3) 



with 5^(x)(a)='| ^'^Q^^jg indicating that A(x) = a is used to 

interpret the voxel xin the test image. Denoting the conditional 
probability of the atias-label selector as PAixMhW"), the Q: 
function reduces to: 



Q(W; W") = £ ,^|, {\ogp(l, W\A)\1, W'''} 
^ P^(,)(fljI(x),W-''''(x)) (4) 

Iog;;(I(x),IF(x)|a) 

The sequence of iterates IF^'',IF*^',..., associated to the alternating 
maximization defined by the iteration 



Step5: Update tiie parcellation W"'''*-W"™ and the optimal 
diffeomorhiphisms (jfj'^ , go to Step 2. 

II — W''''^\f' 
Stop the iteration if either ^ — — <le~ or the 

II W'-if 

number of total iterations is bigger than 30. 

Remarks. 1 . To initiaHze the optimal diffeomorphism ^^"^ in 
Step 1 that is associated to the atias-label pair a, we used the 
optimal diffeomorphism obtained from a two-channel LDDMM 
image mapping with one channel being the FA images and the 
other the MD images, which has been validated in registering DTI 
images [60]. Given the pair of the target /and an atiasJ, we 
compute a diffeomorphic deformation ^between the two vector 
valued images /= [/FA)/MD]and/= [/fa>>^md] such that 
J = Io(j)~^ or [/fAj-^md] = [/fa°i/' './md"'?' ']■ The diffeomorph- 
ism is assumed to be generated as the end point, (p = (j>\, of the flow 
of the smooth time-dependent vector field, V(eK,;e[0,I], via the 
ordinary differential equation ^ = v,(^J'),f6[0,I], where ^^li the 
identity transformation. The optimal diffeomorphic deformation is 
generated by integrating the vector field, which is found to 
minimize the energy: 



V, = arg mm 
1 



p,\\y'dt+ 2ll''FA 

0 <^FA 



1 



/fa°</' ' — /p. 



(9) 



l2 
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where the parameters (TpAand <tmd control the weighting of the 
two contrast-matching terms of smoothness regularization terms. 
In this study, we setffpA = ^MD = 1 ■ To ensure that the solution of 

Eq. (9) lies in the space of difFeomorphisms, the set of time-indexed 
vector fields v, must be sufficiently spatially smooth, requiring Kto 
be a reproducing kernel HUbert space. For computational 
purposes, we use an operator-induced norm on Fsuch that 

11/11 = E 11^/112' and \\fi\\2'= Slfifdx with the differential 
1=1 

operatorL= — aV^^-t-y, where V^^is the Laplacian operator with 
power/>>1.5. In this study, we use/>=l,y=l, and ais selected 
according to the cascading method described in [60,67] 
asO.Ol -0.005-0.002. 

For the initialization of the parcellation label in Step 1 , there are 
multiple choices. In our case, we use the propagation of the labels 
of atlas-label pair ai under the optimal global diffeomorphism 
^{0)...fr'"opW-l. 

2. To incorporate the local optimized dilfeomorphism qfj^'ix) in 
the calculation of Eq. (6), we use mode approximation via 

p{l,W\a)= ^p{l,W\a,(p)Ti{q,\a)d^^p{l,WMa) ■ (10) 

The optimized diffeomorphism in Eq. (10) is obtained as the 
mode, as computed in Eq. (8). 

3. In calculating the terms in Eq. (6), we assume that the prior 
distribution on the adas-label pair^is uni- 
form7c(^ = a)= ,a = a\,...,an. Via Bayes' rule, we have: 



p{I{x),W'"\x)\a4f{x)) 



p(I(x)\ W'^xla^fJ^xM W'"''{x)\a,C(x)) 



(11) 



To calculate;7(I(x)| W'''''{x),a,^f(xj), we define a hierarchical 
model between the image land the underlying diflFeomorphic 
change in coordinates of the atias so that IFspfits land 
Conditioned on W, the joint measurement, I and q>^, is 
independent, giving rise to: 

p(l(x)\ W'''''ix),a,g>fix)) =pilix)\ W'^ixU) . (12) 
Therefore, we have: 

p{I{x),W''"'{x)\a4''^{x)) 

=p(I{x)\ W"'\x),a4f{x))p{W'''''{x)\a4f{x)) (13) 
=p(I{x)\ W''\x),a)p(W"'\x)\a,^f{x)) , 

where p(l(x)\W''''^(x),d)h computed as 

p{l{x)\ W'''''{x),a)=p(lFA{x)JMD(x)Mx)Jy{x)Mx)\ W'''\x),a) 



n p(I„{x)\W°'''(x),a) , 

me{FA,MD,x,y,z} 



(14) 



with I^a(x) : Q^R indicating the FA \^alue at voxel x in the 
target. In calculating each single term in Eq. (14) such as 
P{Ipa(x)\ W^ixXa), we model it as the probability density 
function of a Gaussian Mixture Model (GMM), the parameters 
of which are computed from the atlas-label pair. To be specific, for 
atias-label pair a, we model 



p(lUx)\ W{x) = k)= J2p(Ifa(x)\ Wyx) = k,t)af, (15) 



where denotes the total number of Gaussians in the mixture 
model, p(I§^{x)\W{x) = k,t) represents the probability density 
function of a single Gaussian 



p{I^j,(x)\W(x) = k,t)- 



1 J (I^^(x)-nff \ 



and ^ a"* = 1 are the mixing coefficients for different Gaussians. 
t 

For the parameters of the mixture of Gaussians associated to a 
specific labelA:, 



(17) 



we employ the EM algorithm to derive the maximum-likelihood 
estimators. The termp(/FA(^)| '^'''''(x),a) in Eq. (14) is computed 
according to 



p{lFA{xW'''(x) = k,a) 



M 1 



exp< 



(/fa(x)- 



(18) 



For any given structure, the total numbers of Gaussians, M, for 
the mixtures are pre-defined. We set M = 2 for the structures 
smaller than 1000 mm^ and A/ = 4 for those larger than 1000 
rnm^ . These parameters in the GMM were empirically deter- 
mined. Tlu; GMM is used to c[uantitatively characterize the 
characteristic shape of the histogram of the intensity distribution of 
each contrast in each structure. 

4. To compute the Qjfunction as described in Eq. (7), 
according to Bayes' rule, we have p(l(x),W{x)\a) = 
p(l(x)\W{x),a)p{Wix)\a). The term p{l{x)\W{x),a)vi, computed 
as demonstrated in Eqs. (14) - (18), and p{W{x)\a) is approximated 
via W'>q>~^ under trUinear interpolation. 

5. Given our splitting assumption in Eq. (12), Eq. (8) is 
equivalent to q)^^^ = a.rgTaaxp(lV"™\a,ip)n{ip\a). Considering 

computational efficiency, we use measures of the distance between 
the parcellation of the target and the difFeomorphically deformed 
results of the atias parcellations, analogous to LDDMM for image 
matching and surface matching. To be specific, we use the Dice 
overlap measurement between IF"™ and IF"o^~' to approximate 
the te.rmp(W"™\a,(p). The optimal local diffeomorphism is 
assumed to come from a composition of the optimal global 
diffeomorphism and an optimal local 12-parameter affine trans- 
formation. Namely, (pjix) = qf^^°a{x), where qf-^^ is the optimal 
global diffeomorphism computed in Step 1 and ais the optimal 
local 12-parameter affine transformation that maximizes 
;)(IF"™|a,^*'"oa)7t(^(Ooa|a), where /)(IF"''»'|a,^<,0'oa)is quantified 
as the Dice overlap between ir"™and W''°{q>^^^°ay^ . Note that 
the optimal local affine transformation is obtained on a structure- 
by-structure basis. Therefore, for a single atlas-label pair a, the 
optimized local diffeomorphisms (p"™{x) should be identical for 
voxels in the same structure. Given atias-label pair a, the prior 
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Figure 1. Demonstration of the unique anatomical features revealed by multi-contrast images generated in DTI and GMM. 

Histograms of the five contrasts, FA, IVID, EV-x, EV-y, and EV-z, of five adjacent structures are shown, including two white matter structures (the ALIC 
and the PLIC), two gray matter structures (the caudate and the thalamus), and the ventricle. In each subplot, blue indicates the histogram of the 
corresponding contrast within that specific structure, green represents the probability density of each single Gaussian, and red shows the weighted 
sum of all Gaussians. Abbreviations are: ALIC: Anterior limb of internal capsule and PLIC: Posterior limb of internal capsule. 
doi:1 0.1 371/journal.pone.0096985.g001 



distribution of the transformation n{(p\a) is estimated as tlie 
multiplication of two terms n((p\a) = n(q/^''\a)n(oi\a,q)^^^), wliere 
7i(^*J''|fl)is estimated by the one over the metric distance [68] in 
diffeomorphism space given by tlie exponential of the geodesic 
length, computed from the two-channel LDDMM mapping. The 
prior on the 12-parameter affme transformation7i(o(|a,^Jj''') is 
modeled as a multivariate Gaussian, n(oi\a,^'^^^) = N{Mai,Cai), 
similar to the strategy adopted in [69]. In our approach, we 
use M^ = [l,0,0,0,l,0,0,0,l,0,0,0]^andQ = diag([le-Me-2 
le^^,le^,le^,le^]). We assume that all the parameters are 
mutually independent, and thus the covariance matrix is diagonal. 
Since the first 9 parameters in arepresent the affine matrix, their 
variances should be small, for which we assign 0.01. The last 3 
parameters represent the translation in the :\:,3^,zdirections, and 
therefore their variances should be big, for which we use 100. 

Since the optimization of the local diffeomorphisms is based on 
the overlap between the parcellation W^'^oi the target and the 
diffeomorphicaUy deformed results of the adas parcellations and 
the optimized diffeomorphisms ^„{x) are identical for voxels in the 
same structure, the term p(lV''''^(x)\a,q>°'''(x))m Eq. (11) is 



approximated as being proportional to the overlap distance 
between W"'''{x) and W''o(p<>J''{- \){x). Again, for adas-label pair 
a, this quantity is identical for voxels in the same structure. 

To sum up, the MAP estimation problem is solved in an EM 
approach. We iterate between fixing the local optimal diffeo- 
morphism for each label in each atlas-label pair and obtaining the 
maximizing parcellation of the target, and then locally optimizing 
the diffeomorphisms associated to each label in each atlas-label 
pair given the frxed parcellation. 

Image quantification 

After the brain had been parceUated to the 159 structures, the 
peripheral ROIs were further decomposed to the CSF, cortex, and 
peripheral white mater using MD (threshold value = 0.0015 to 
separate the CSF and the tissue) and FA (threshold value = 0.2 to 
separate the cortex and the white matter) (Faria et al., 2010; Oishi 
et al., 2009). The CSF regions were excluded from the analysis. 
There were 50 peripheral ROIs and thus the final number of 
ROIs was 193. The volumes of these ROIs were obtained by 
counting the number of voxels. ROI-specific FA and MD values 
were measured by averaging the values of all voxels within the 
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Figure 2. A comparison of parceliating using a single-contrast image and multi-contrast images, in terms of overlap accuracy. Tine 
mean Dice overlaps and the standard deviations of the eighteen ROIs obtained from automated parcellations based on five contrasts (red), the single 
FA contrast (green), the single MD contrast (blue), the vector x y z contrasts (yellow), as well as the inter-rate (patterned). The mean values are 
calculated across fourteen different subjects. Star marks indicate significant difference among the four sets of Dice results by ANOVA (p<<0.05). 
Abbreviations are: GCC - genu of corpus callosum; BCC - body of corpus callosum; Caud - caudate; Put - putamen; ALIC - anterior limb of internal 
capsule; PLIC - posterior limb of internal capsule; CG - cingulate gyrus; IVICP - middle cerebellar peduncle; SLF - superior longitudinal fasciculus; 
CST - corticospinal tract. 
doi:1 0.1 371 /journal.pone.0096985.g002 
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Figure 3. A correlation comparison of the automated caudate parcellation from using a single-contrast image and multi-contrast 
images. A plot of the correlation between the automated and the manual measurements of the size of the caudate in both hemispheres in square 
millimeters. Results from the four automated parcellation methods are compared: 5-contrast (red), FA-only (green), MD-only (blue), and EV-only 
(yellow). 

doi:1 0.1 371 /journal.pone.0096985.g003 
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Figure 4. A comparison of the CST correlation obtained from using a single-contrast image and multi-contrast images. A correlation 
plot between the automated and manual measurements of the sizes of left and right corticospinal tracts (CST). Results from the four automated 
parcellation methods are compared: 5-contrast (red), FA-only (green), IVID-only (blue), and EV-only (yellow). 
doi:1 0.1 371 /journal.pone.0096985.g004 



ROIs. The parcellation criteria used in this paper foUowed our 
previous pubhcations [70,71], in which the cortex and the white 
matter definitions followed ICBM-LPBA40 [24] and probabilistic 
white matter adas [71], respectively. 

Manual delineation for accuracy measurements 

Eighteen structures (sixteen white matter structures and two 
deep gray matter structures) were manually delineated on the pre- 
selected 2D slices of fourteen subjects from three groups - four 
from the normal group, five from the mild abnormal group, and 
five from the severe abnormal group. The manual delineation was 
performed by incorporating information from MD, FA, and color- 
coded eigenvector maps. To investigate the intra- and inter-rater 
variability of the manual delineations, two raters pC.T. andJ.H.) 
performed the manual delineations twice with more than 3-week 
intervals. To quantitatively evaluate the parcellation accuracy of 
our algorithm, we used four measurements: 

1. Dice overlap coefficients 

We calculated the Dice overlap coefficients between the 
manually delineated 2D ROI and the corresponding ROI in the 
automated parcellations. The Dice overlap coefficient is calculated 
2TP 

as: D = , where TP is the area of the region that 

2TP + FP + FN ^ 

belongs to both the automated ROI and the manual ROI, FP is 
the area of the region that belongs to the automated ROI but not 
the manual, and FN is the area of the region that belongs to the 
manual ROI but not the automated. 

2. The correlation between the size of the manually delineated 
ROI and that of the automated ROI. 



3. The correlation between the mean FA value of the manual 
ROI and the mean FA of the automated ROI. 

4. The correlation between the mean MD value of the manual 
ROI and the mean MD of the automated ROI. 

To evaluate the improvement in parcellation accuracy given by 
the multi-contrast approach, we compared the parcellations from 
the 5-contrast multi-atlas approach (FA, MD, vector elements x, y, 
and z, combined), with those obtained from the same multi-adas 
but with only a single contrast — FA-only, MD-only, and a three- 
contrast approach — EV-only (x, y, z combined). These four 
methods vary from each other only in the computation of Eq. (14). 
To compare the four methods statistically, for each structure, we 
performed a one-way ANOVA to examine significant difference 
among the Dice results obtained from the four approaches. For 
statistical correlation analysis, we used WUliam's modification of 
Hotelling's test [72]. 

For the scan-rescan reproducibility test, we investigated 
the volume difference between the automated parcellations 
of the same structure from the two scans for the 
same subject. The volume difference is computed as: 
|vo/(ROIi)-vo/(ROl2)| 



VD-- 



-, where vo/(ROIi)denotes the 



(vo/(ROIi) + vo/(ROl2))/2' 
volume of a specific ROI for scan 1 and vo/(ROl2)denotes the 
volume of the same ROI for the second scan of the same subject. 
In addition, we examined the difference between the mean FA 
value of the automated parcellation of each single structure for die 
first scan and that of the automated parcellation for the second 
scan, as well as the difference between the mean MD values. 
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Results 

Figure 1 demonstrates the concept of the multi-contrast image 
parcellation using five different contrasts obtained from DTI as 
well as the concept of characterizing the intensity distribution of 
each contrast using a GMM. The five selected structures are 
spatially adjacent to each other and need to be accurately 
demarcated based on their contrast features. Each of the five 
contrast rows cannot uniquely differentiate all the five structures, 
but each column (structure) has a unique contrast signature by 
combining these five contrasts. For example, the distinction 
between the tissue and the ventricle is most clear on the MD 
image, while the distinction between the caudate and the anterior 
limb of internal capsule (ALIO) is most clear on the FA image. 
Likewise, the difference between ALIO and the posterior limb of 
internal capsule (PLIC) is largest in the eigenvector image. The 
GMM quantitatively characterizes the contrast features delineated 
by these multi-channel histograms. 

Figure 2 shows results of Dice measurements, reporting spatial 
agreement with manual delineation. The 5-contrast approach is 
compared with FA-only, MD-only, and EV-only approaches. 
Because of the unique contrast signature of each structure, the best 
contrast that can accurately define it varies. For example, to define 
the contricospinal tract (CST), the EV provides the best accuracy, 
but it provides poor results to define the putamen, which is best 
defined by FA or MD. The 5-contrast approach performs well for 
all structures. According to the results from the one-way ANOVA, 
we found statistical differences among the 4 approaches in 1 1 
structures, in which the 5-contrast approach consistently achieved 
one of the best results. These structures include: the caudate, the 
putamen, the cingulate gyrus, the middle cerebellar peduncle, and 
the corticospinal tract in both hemispheres. The absolute Dice 
level was 0.8-0.9. Note that some structures are difficult to define 
even manually with perfect reproducibility. A good example is the 
superior longitudinal fasciculus (SLF), which has a vague structural 
boundary and the inter-rater variability is large (Dice = 0.8+/ — 
0.259). Because the manual defmition is used as the gold standard, 
the spatial matching cannot be better than the inter-rater spatial 
matching (automated methods cannot be more accurate than 
manual delineation by definition). 

The correlation coefficients between the sizes of the manual and 
the automated parcellations obtained from the four approaches for 
all the eighteen ROIs are listed in Table 3. For some structures, 
the ROI sizes from aU the four automated approaches are highly 
correlated with the ROI sizes of the manual delineations. 
However, structures such as the caudate, the corticospinal tract 
(CST), and the cingulate gyrus (CG), the performance varies from 
approach to approach. In Figure 3 and Figure 4, we show 
examples of the correlation plot between the automated and 
manual approaches for a gray matter (caudate) and a white matter 
(CST) structures. 

Figure 5 shows actual parcellation results of the CST in the 
brainstem of subjects with different degrees of abnormalities, 
which demonstrates how the integration of five contrast informa- 
tion can accurately delineate the sizes. In this example, the fiber- 
orientation information in the EV contrast is necessary to 
accurately reflect the small CST sizes in Case #3. Namely, the 
CST has a characteristic Z-orientation (blue) fiber orientation, 
which can uniquely differentiate the CST from the surrounding 
high-FA white matter structures. The integration of the EV 
information provides strong constraints for the parcellation, 
specifically defining the high-FA regions with a strong orientation 
alignment along the Z axis. The FA-, MD-, and EV-only 
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Figure 5. Examples of CST parcellations from single- and multi-contrast approaches. Demonstration of the parcellation accuracy of the 
CST in three representative cases with different degrees of anatomical abnormalities. Results from five different approaches are compared. 
doi:10.1371/journal.pone.0096985.g005 



approaches extracted the CST accurately for Case 7^1 and ^2, 
but grossly overestimated the CST size for Case #3. 

Based on the comparison results shown in Table 3, significant 
improvement of the correlation between the size of the automated 
parcellations and that of the manual delineations was achieved by 
the 5-contrast approach, compared with the other single contrast 
approaches. Likewise, we performed the manual-auto correlation 
analyses of the mean FA and MD values within each single ROI. 
As shown in Table 4 and 5, again, the 5-contrast approach is 
consistently superior to the other three approaches in terms of 
either FA or MD correlation. 

Figure 6 demonstrates parcellation results for three patients with 
different degrees of abnormalities. A high level of parcellation 
accuracy is visually appreciable for the wide variety of anatomical 
states. 

According to our test-retest experiments, the reproducibility 
results were 4.7%, 2.19%, and 2% for the volume, FA, and MD, 
respectively, averaged over all 193 structures. If we remove 26 
small gray matter structures (<1000 mm'), the reproducibility 
improves to 3.73%, 1.91%, and 1.79% respectively. These small 
gray matter structures had poor test-retest reproducibility because 
they lack clear contrasts in DTI and they are usually not the 
targets of DTI measurements. 

Figure 7 shows the maps of cross-subject variability in the 
volume, FA, and MD measurements. The cross-subject variability 
is computed for each label of interest. It is quantified as the ratio of 
the standard deviation to the mean value across the sixteen 
subjects. A large amount of morphological variability was found in 
the ventricle volumes, while the standard deviations of the volumes 
of white matter structures are in 10-20% range. The standard 
deviations of the FA and MD were noticeably lower and most 
areas were below 10%. The table in Appendix shows a 
comprehensive report of the test-retest reproducibility and the 
cross-subject variability of all the 193 defined structures. These 
values should provide useful information for power calculation in 
future study designs. 

Discussion 

In this study, we developed and tested an automated image 
parcellation method based on a multi-contrast multi-atias likeli- 
hood-fusion algorithm. DTI can generate multiple quantitative 



maps with markedly different qualities of anatomical contrasts. 
The mean diflfusivity contrast provides clear distinction between 
the tissue (generally within the range of 0.6-0.9 cmVs) and the 
CSF (approximately 3.0-3.5 cm^/s), providing a strong con- 
straints to define the ventricles and the brain surface. The FA 
contrast provides sharp distinction between the gray (typically 
FA<0.15-0.25) and white matter structures (FA>0. 15-0.25). The 
eigenvector (EV) can differentiate intra-white matter structures 
based on their characteristic orientations. In this work, we used the 
absolute values of the three components, EV-x, EV-y, and EV-z. 
While this approach solves the difficulties associated with the sign 
of the eigenvectors, some orientation information degenerates 
[73]. This is obviously a simplified approach and there is room for 
improvement. The mixture of these three types of information 
could also invite noise. For example, in the low-FA gray matter 
structures, the fiber orientation information may be random and 
should not receive significant weighting. This type of weighting is 
naturally achieved by incorporating the variability information 
about the voxel values within a single parceUated structure; for 
example, the EV information of the thalamus in Figure 1 shows 
almost equal values for the X, Y, and Z channels with high intra- 
structure variability. In Eq. (18), we use mixtures of Gaussians to 
model the intensity distribution within a single structure. If the 
intensity within the structure is homogeneous, the algorithm 
automatically assigns weight 1 to a single Gaussian. If there is high 
intra-structure variability, multiple Gaussians will be used to 
model the intensity, with each Gaussian being given a small 
weight. In this way, it effectively reduces the contribution of this 
contrast information in computing the quantity 
p{I„,(x)\ fF,fii),me{FA,MD,EV - x,EV - y,EV - z} in Eq. (14). 
By incorporating the consistent anatomical signatures into the 
parcellation criteria, we aim to achieve robust parcellation. In the 
past, multi-contrast image registration approaches have been 
postulated including ones for DTI data [60,74-79]. These tensor- 
or vector-based registration also indicated improved registration 
accuracy [56-59] . The proposed method can be considered as an 
extension of these previous works by incorporating them into a 
multi-atlas framework. 

The improved accuracy, with respect to a single-contrast 
approach, is shown in Figure 2 using Dice measurements. While 
the performance of single-contrast approaches varies depending 
on the structure, the five-contrast approach coiisistendy achieved 
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Figure 6. Examples of whole brain parcellations. Results of the whole brain parcellations into 159 structures in three representative cases with 
large anatomical variability. The parcellation results are superimposed on color (upper row) and MD (bottom row) images. 
doi:1 0.1 371 /journal.pone.0096985.g006 

the highest level of accuracy. As mentioned in the Results section, good as a human rater. Careful observation of this figure reveals 
the accuracy of the automated method cannot be higher than the that, for many core white matter structures with distinctive and 
reproducibility level of the manual delineation. In this sense, the vmiform tract orientations, the eigenvector contrast alone can 
results in Figure 2 indicate that the five-channel approach is as provide a similar level of accuracy as the five-contrast method. 




Figure 7. Depiction of the cross-subject variability in different whole brain structures of interest. Demonstration of the cross-subject 
variability (std/mean) within the 16 healthy subjects for each of the 193 anatomical regions. The population variability in terms of ROI volume (top 
row), the mean MD of each ROI (middle row), and the mean FA of each ROI (bottom row) are shown. The results are presented for three difference 
axial slices using a colormap with the color scale ranging from 0 to 0.5. 
doi:1 0.1 371/journal.pone.0096985.g007 
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suggesting that the parcellation was mainly driven by the 
eigenvector contrast. However, if a structure is not characterized 
by a uniform fiber orientation due to curvature within a segment, 
such as GCC and BCC, the eigenvector may not always be the 
most reliable contrast. 

While the Dice measurements (Figure 2) provide important 
information about the accuracy level of the automated parcella- 
tion, it is probable that the correlation results reported in Table 3 
are more important for actual image-based studies. Anatomiced 
delineations of brain structures depend on anatomical definitions. 
It is reasonable that there is consistent difference in boundaries of 
defined structures between two different approaches. In this sense, 
low Dice vEilues do not necessarily mean that the automated results 
are not useful. The high correlation between the manual and 
automated methods indicates that they have similar powers to 
differentiate different anatomical states, which is ultimately the 
goal of quantitative analyses. In this respect, the high correlation 
between the 5-contrast and manual approaches is an encouraging 
result. 

One limitation of our accuracy evaluation is that the measured 
structures were limited to the core brain structures, which can be 
reproducibly defined manually. This is inevitable because the 

manual delineation results were used as the gold standard. As 
reported in the Results section, one of the core white matter 
structures, the SLF, suffered from low inter-rater reproducibility 
due to its complex shape. Reproducible definitions of peripheral 
white matter regions by manual tracing would be prohibitively 
difficult and, due to the absence of the gold standard, the accuracy 
measurements of the proposed automated segmentation were 
challenging. In this study, we are therefore limited to reporting 
test-retest reproducibility and population variability measure- 
ments, which could be important resources for power analyses of 
the proposed method. 

The test-retest reproducibility showed less than 5% variability 
for the volume measurement for most of the defined structures (see 
Appendix table). The test-retest reproducibility measures of FA and 
MD indicated higher reproducibility (less than 3%). Th(" 
anatomical variability for the 193 measured structures reported 
in Figure 7 should provide information for power analysis to 
design population studies. 

In this study, we reported the accuracy level of the multi- 
contrast multi-atlas apjjroacli for a wide range' of anatomical 
phenot)'pes. The performance' of this technology relies heavily on 
the availability of atlases with (;:onsistent parcellation criteria. 
Creation of such atlases is a time consuming task, which is one of 
the limitations of this approach. The accuracy of the parcellation 
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